03_models

loading packages

suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(sjPlot))
suppressPackageStartupMessages(library(lme4))
suppressPackageStartupMessages(library(performance))
suppressPackageStartupMessages(library(nortest))

options(scipen = 999)

importing dataset

dat_model <- read_csv("data/data_model.csv")

baseline model

trait rumination

# building model
trait_model <-
  lmer(
    NA_sum ~ age + gender + beep + day +
    event_unpleasantness_centered * trait_rumination +
    (event_unpleasantness_centered | id),
      data = dat_model
    )


# model outputs
summary(trait_model)

Linear mixed model fit by REML ['lmerMod']
Formula: NA_sum ~ age + gender + beep + day + event_unpleasantness_centered *  
    trait_rumination + (event_unpleasantness_centered | id)
   Data: dat_model

REML criterion at convergence: 33508.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.5777 -0.4677 -0.0848  0.2800  7.6410 

Random effects:
 Groups   Name                          Variance Std.Dev. Corr
 id       (Intercept)                   0.77855  0.8824       
          event_unpleasantness_centered 0.03429  0.1852   0.59
 Residual                               0.56096  0.7490       
Number of obs: 14265, groups:  id, 247

Fixed effects:
                                                 Estimate Std. Error t value
(Intercept)                                     0.2160764  0.1728315   1.250
age                                            -0.0047372  0.0041325  -1.146
genderMale                                     -0.1152025  0.1390809  -0.828
genderOther                                     0.0497428  0.7339194   0.068
beep                                            0.0051633  0.0029599   1.744
day                                            -0.0047750  0.0008425  -5.668
event_unpleasantness_centered                   0.3053063  0.0143028  21.346
trait_rumination                                0.3551413  0.0664066   5.348
event_unpleasantness_centered:trait_rumination  0.0779311  0.0150512   5.178

Correlation of Fixed Effects:
            (Intr) age    gndrMl gndrOt beep   day    evnt__ trt_rm
age         -0.925                                                 
genderMale  -0.190  0.040                                          
genderOther -0.071  0.047  0.059                                   
beep        -0.069 -0.001 -0.002  0.000                            
day         -0.041 -0.002  0.002 -0.003  0.008                     
evnt_nplsn_  0.159  0.011 -0.036  0.013  0.035  0.011              
trait_rmntn -0.308  0.281  0.289  0.082 -0.006 -0.005 -0.014       
evnt_npl_:_  0.010 -0.018  0.030 -0.015  0.003 -0.003 -0.078  0.457

tab_model(trait_model, show.std = T)


# assumption checks
check_collinearity(trait_model)

# Check for Multicollinearity

Low Correlation

                                           Term  VIF       VIF 95% CI Increased SE Tolerance Tolerance 95% CI
                                            age 1.12 [1.11,     1.15]         1.06      0.89     [0.87, 0.90]
                                         gender 1.12 [1.10,     1.14]         1.06      0.89     [0.87, 0.91]
                                           beep 1.00 [1.00,   138.99]         1.00      1.00     [0.01, 1.00]
                                            day 1.00 [1.00,      Inf]         1.00      1.00     [0.00, 1.00]
                  event_unpleasantness_centered 1.01 [1.00,     1.05]         1.01      0.99     [0.95, 1.00]
                               trait_rumination 1.58 [1.54,     1.61]         1.26      0.63     [0.62, 0.65]
 event_unpleasantness_centered:trait_rumination 1.34 [1.31,     1.36]         1.16      0.75     [0.73, 0.76]

plot_model(trait_model, type = "diag")

[[1]]

[[2]]
[[2]]$id


[[3]]

[[4]]

#check_model(trait_model)


# visualizing predicted values
plot_model(trait_model)

# goodness of fit indicators
AIC_trait <- AIC(trait_model)
BIC_trait <- BIC(trait_model)

coefplot for trait rumination

coefplot_trait <- 
  plot_model(trait_model,
             type = "std",
             rm.terms = c("age", "gender [Male]", "gender [Other]", "beep", "day", "event_unpleasantness_centered"),
             ci.lvl = 0.95,
             std.est = T,
             title = "",
             show.values = TRUE,
             colors = "black",
             axis.labels = c("Trait rumination * Perceived stress", "Trait rumination"),
             axis.title = "Standardized β coefficients"
             ) +
             ylim(0, 0.6) +
             theme_minimal() +
             theme(axis.line = element_line(size = .3),
                   axis.text = element_text(size = 10, color = "black"),
                   axis.title.x = element_text(size = 12, color = "black")
                   )

Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

coefplot_trait

mean state rumination

# building model
mean_state_model <-
  lmer(
    NA_sum ~ age + gender + beep + day +
    event_unpleasantness_centered * mean_state_rumi_grand_centered +
    (event_unpleasantness_centered | id),
    data = dat_model
    )

# model outputs
summary(mean_state_model)

Linear mixed model fit by REML ['lmerMod']
Formula: NA_sum ~ age + gender + beep + day + event_unpleasantness_centered *  
    mean_state_rumi_grand_centered + (event_unpleasantness_centered |      id)
   Data: dat_model

REML criterion at convergence: 33378.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.5871 -0.4680 -0.0864  0.2809  7.6518 

Random effects:
 Groups   Name                          Variance Std.Dev. Corr
 id       (Intercept)                   0.42164  0.6493       
          event_unpleasantness_centered 0.03091  0.1758   0.51
 Residual                               0.56109  0.7491       
Number of obs: 14265, groups:  id, 247

Fixed effects:
                                                               Estimate Std. Error t value
(Intercept)                                                   0.1402751  0.1291530   1.086
age                                                          -0.0032922  0.0031130  -1.058
genderMale                                                   -0.0501490  0.1062399  -0.472
genderOther                                                   0.1741375  0.5730614   0.304
beep                                                          0.0049676  0.0029591   1.679
day                                                          -0.0048631  0.0008409  -5.783
event_unpleasantness_centered                                 0.3068529  0.0138163  22.209
mean_state_rumi_grand_centered                                0.9105524  0.0598964  15.202
event_unpleasantness_centered:mean_state_rumi_grand_centered  0.1227505  0.0188286   6.519

Correlation of Fixed Effects:
            (Intr) age    gndrMl gndrOt beep   day    evnt__ mn____
age         -0.922                                                 
genderMale  -0.129 -0.027                                          
genderOther -0.053  0.028  0.045                                   
beep        -0.095 -0.001 -0.002  0.001                            
day         -0.055 -0.004  0.005 -0.004  0.007                     
evnt_nplsn_  0.139  0.005 -0.033  0.010  0.036  0.011              
mn_stt_rm__ -0.182  0.155  0.200  0.061 -0.009 -0.001 -0.015       
evnt__:____  0.016 -0.021  0.001 -0.012  0.004 -0.002 -0.063  0.387

tab_model(mean_state_model, show.std = T)


# assumption checks
check_collinearity(mean_state_model)

# Check for Multicollinearity

Low Correlation

                                                         Term  VIF       VIF 95% CI Increased SE Tolerance
                                                          age 1.04 [1.02,     1.06]         1.02      0.96
                                                       gender 1.06 [1.04,     1.08]         1.03      0.94
                                                         beep 1.00 [1.00,    64.69]         1.00      1.00
                                                          day 1.00 [1.00,      Inf]         1.00      1.00
                                event_unpleasantness_centered 1.01 [1.00,     1.08]         1.00      0.99
                               mean_state_rumi_grand_centered 1.28 [1.26,     1.31]         1.13      0.78
 event_unpleasantness_centered:mean_state_rumi_grand_centered 1.20 [1.18,     1.23]         1.10      0.83
 Tolerance 95% CI
     [0.94, 0.98]
     [0.92, 0.96]
     [0.02, 1.00]
     [0.00, 1.00]
     [0.93, 1.00]
     [0.76, 0.79]
     [0.82, 0.85]

plot_model(mean_state_model, type = "diag")

[[1]]

[[2]]
[[2]]$id


[[3]]

[[4]]

#check_model(mean_state_model)


# visualizing predicted values
plot_model(mean_state_model)

# goodness of fit indicators
AIC_mean_state <- AIC(mean_state_model)
BIC_mean_state <- BIC(mean_state_model)

coefplot for mean state rumination

coefplot_mean_state <- 
  plot_model(mean_state_model,
             type = "std",
             rm.terms = c("age", "gender [Male]", "gender [Other]", "beep", "day", "event_unpleasantness_centered"),
             ci.lvl = 0.95,
             std.est = T,
             title = "",
             show.values = TRUE,
             colors = "black",
             axis.labels = c("Mean state rumination * Perceived stress", "Mean state rumination"),
             axis.title = "Standardized β coefficients"
             ) +
             ylim(0, 0.6) +
             theme_minimal() +
             theme(axis.line = element_line(size = .3),
                   axis.text = element_text(size = 10, color = "black"),
                   axis.title.x = element_text(size = 12, color = "black")
                   )

Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

coefplot_mean_state

state rumination

# building model
state_model <-
  lmer(
    NA_sum ~ age + gender + beep + day +
    event_unpleasantness_centered * rumination_centered +
    (event_unpleasantness_centered | id), data = dat_model
    )


# model outputs
summary(state_model)

Linear mixed model fit by REML ['lmerMod']
Formula: NA_sum ~ age + gender + beep + day + event_unpleasantness_centered *  
    rumination_centered + (event_unpleasantness_centered | id)
   Data: dat_model

REML criterion at convergence: 31593.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.5625 -0.4499 -0.0771  0.2846  8.2712 

Random effects:
 Groups   Name                          Variance Std.Dev. Corr
 id       (Intercept)                   0.86414  0.9296       
          event_unpleasantness_centered 0.02793  0.1671   0.66
 Residual                               0.48929  0.6995       
Number of obs: 14265, groups:  id, 247

Fixed effects:
                                                    Estimate Std. Error t value
(Intercept)                                        0.4238103  0.1637515   2.588
age                                               -0.0090868  0.0039334  -2.310
genderMale                                        -0.2758397  0.1338196  -2.061
genderOther                                       -0.1560593  0.7273910  -0.215
beep                                               0.0009888  0.0027650   0.358
day                                               -0.0045018  0.0007871  -5.720
event_unpleasantness_centered                      0.2545536  0.0129058  19.724
rumination_centered                                0.2651076  0.0066691  39.751
event_unpleasantness_centered:rumination_centered  0.0386275  0.0038579  10.013

Correlation of Fixed Effects:
            (Intr) age    gndrMl gndrOt beep   day    evnt__ rmntn_
age         -0.912                                                 
genderMale  -0.094 -0.066                                          
genderOther -0.039  0.016  0.031                                   
beep        -0.072  0.002  0.001  0.001                            
day         -0.043  0.001  0.002 -0.002  0.008                     
evnt_nplsn_  0.206  0.003 -0.050  0.012  0.040  0.011              
rmntn_cntrd  0.015 -0.009 -0.014 -0.005 -0.029 -0.005 -0.076       
evnt_npl_:_ -0.009  0.001  0.006 -0.001 -0.005  0.022 -0.031 -0.299

tab_model(state_model, show.std = T)


# assumption checks
check_collinearity(state_model)

# Check for Multicollinearity

Low Correlation

                                              Term  VIF       VIF 95% CI Increased SE Tolerance Tolerance 95% CI
                                               age 1.00 [1.00,     1.15]         1.00      1.00     [0.87, 1.00]
                                            gender 1.01 [1.00,     1.07]         1.00      0.99     [0.94, 1.00]
                                              beep 1.00 [1.00,     2.67]         1.00      1.00     [0.37, 1.00]
                                               day 1.00 [1.00, 1.87e+07]         1.00      1.00     [0.00, 1.00]
                     event_unpleasantness_centered 1.01 [1.00,     1.05]         1.01      0.99     [0.96, 1.00]
                               rumination_centered 1.11 [1.09,     1.13]         1.05      0.90     [0.88, 0.92]
 event_unpleasantness_centered:rumination_centered 1.10 [1.08,     1.12]         1.05      0.91     [0.89, 0.92]

plot_model(state_model, type = "diag")

[[1]]

[[2]]
[[2]]$id


[[3]]

[[4]]

#check_model(state_model)


# visualizing predicted values
plot_model(state_model)

# goodness of fit indicators
AIC_state <- AIC(state_model)
BIC_state <- BIC(state_model)



# outputs for trait, mean state and state models
tab_model(trait_model, mean_state_model, state_model, show.std = T)

coefplot for state rumination

coefplot_state <- 
  plot_model(state_model,
             type = "std",
             rm.terms = c("age", "gender [Male]", "gender [Other]", "beep", "day", "event_unpleasantness_centered"),
             ci.lvl = 0.95,
             std.est = T,
             title = "",
             show.values = TRUE,
             colors = "black",
             axis.labels = c("State rumination * Perceived stress", "State rumination"),
             axis.title = "Standardized β coefficients"
             ) +
             ylim(0, 0.6) +
             theme_minimal() +
             theme(axis.line = element_line(size = .3),
                   axis.text = element_text(size = 10, color = "black"),
                   axis.title.x = element_text(size = 12, color = "black")
                   )

Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

coefplot_state

full model

check_model(full_model)

Package `see` required for model diagnostic plots.

coefplots for full model

# plot for interactions (colored afterwards in photoshop)
coefplot_interactions <- 
  plot_model(full_model,
             type = "std",
             rm.terms = c("rumination_centered", "trait_rumination", "mean_state_rumi_grand_centered", "age", "gender [Male]", "gender [Other]",   "beep", "day", "event_unpleasantness_centered"),
             ci.lvl = 0.95,
             std.est = T,
             title = "",
             show.values = TRUE,
             colors = "black",
             axis.labels = c("Trait rumination * Perceived stress", "Mean state rumination * Perceived stress", "State rumination * Perceived stress"),
             axis.title = "Standardized β coefficients"
             ) +
             ylim(0, 0.1) +
             theme_minimal() +
             theme(axis.text = element_text(size = 10, color = "black"),
                   axis.title.x = element_text(size = 12, color = "black"),
                   axis.line = element_line(size = .3)
                   )

Warning: Model failed to converge with max|grad| = 0.00403381 (tol = 0.002, component 1)Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

coefplot_interactions

# plot for main effects (colored afterwards in photoshop)
coefplot_maineffects <- 
  plot_model(full_model,
             type = "std",
             terms = c("rumination_centered", "mean_state_rumi_grand_centered", "trait_rumination"),
             ci.lvl = 0.95,
             std.est = TRUE,
             title = "",
             show.values = TRUE,
             colors = "black",
             axis.labels = c("Trait rumination", "Mean state rumination", "State rumination"),
             axis.title = "Standardized β coefficients"
             ) +
             ylim(0, 1) +
             theme_minimal() +
             theme(axis.text = element_text(size = 10, color = "black"),
                   axis.title.x = element_text(size = 12, color = "black"),
                   axis.line = element_line(size = .3)
                   )

Warning: Model failed to converge with max|grad| = 0.00403381 (tol = 0.002, component 1)Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

coefplot_maineffects

saving full model for bootstrapping

saveRDS(full_model, "saved_models/full_model.rds")